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ABSTRACT 


An analytical framework was constructed to interchange ice resistance components from 
existing ice resistance calculation methods. Within this framework, the Lindqvist 
analytical method bending ice resistance component was substituted with a value 
obtained from a finite element static model of the ice sheet. A parametric study of this 
substitution was performed with error correction to within 12% over the entire range of 
parameter variation. The finite element ice sheet model was then damaged; and then the 
average ice bending resistance was obtained, substituted into the Lindqvist analytical 
method, and quantified as a change in the total ice resistance. Therefore, a hybrid ice 
resistance model was developed that accounts for the effect of damage to the bending 
resistance component, and enables further study of unconventional icebreaking methods. 
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I. INTRODUCTION 

A. ICEBREAKING HISTORY 

Prior to the 20th century, the first true icebreaker, Yermak, was commissioned 
and tasked with exploring the polar regions of the globe. The vessel had a displacement 
of 10,000 tons, a 10,000 hp propulsion system consisting of three screws aft with one 
screw forward, thickened hull plating, double hull construction, a ballast system to 
quickly affect roll/trim, and a rounded hull form [1], The summation of these design 
characteristics ensured this to be the first vessel to reliably transit multi-year sea ice, 
which it did from 1899 to 1963. The success of this ship ultimately defined the 
specifications for the lACS Polar Class of vessels, from which all modem icebreakers are 
derived [2]. As a result of this tried and tme formula for immensely powerful and 
specialized ice navigation platforms, the overall ship design and mechanical process to 
break ice is essentially the same today as it was back then. Despite this similitude, the 
icebreaker underwent several design iterations, and efficiency is considerably improved 
over 1899 as the concept was refined and matured. 

Substantial research was also conducted to develop systems to damage ice ahead 
of the icebreaker such that the ice failed with less effort from the ship. This concept was 
successfully used to aid early expeditions to the Antarctic, by fracturing the ice with 
explosives. Less violent means of damaging ice, such as water jet cutting mechanisms, 
were conceived to provide the benefit of the explosives with the practicality of a 
shipboard system. The Yermak served as the test bed for the first water jet cutters aboard 
icebreakers sometime after the idea was introduced in 1935 [3]. The water jets at that 
time lacked the power density and peak pressures to provide any benefit to the 
icebreaker. At the time the Yermak was retired in 1963, follow on research attempted to 
physically model the effectiveness of an icebreaker with ice cutting systems [3]. The 
significant disparity between the failure characteristics of model and full scale ice 
resulted in criticism that the claims of this report were too optimistic in favor of ice 
cutting systems, and more specifically, water jet cutters [4]. Nonetheless, the greatest 
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contribution of this work was the desire to use modeling methods to craft new strategies 
to damage the iee and allow the passage of a ship. 

From literature review, two themes following the 1963 study beeame apparent. 
First, the efficieney of the iee eutting systems received detailed study [4]-[8].Seeond, 
mueh researeh was done to adequately define the material properties of iee [9], [10], and 
apply this to reasonably model the ieebreaking proeess [11], [12]. With minor exeeption 
[3], the researeh tended to ignore the applieation of these ieebreaking systems in terms of 
a grander ieebreaking strategy. In other words, the speeifie damage pattern that eauses iee 
to most readily fail in the presence of the ieebreaker received little attention. Just as the 
summation of its parts cemented the Yermak’s title as the first true icebreaker, the 
effectiveness of the iee breaking / ice cutting system relies on the synergy between the 
available tools and their artful manipulation. 

B. STUDY OBJECTIVE 

This study eombines the advantages of existing iee resistanee ealeulation methods 
into a single hybrid method, such that the steady state veloeity of an ieebreaker is 
eorrelated to the available thrust; and inelude eases in whieh the iee sheet possesses some 
artifieial damage. In addition, this study intends to demonstrate the utility of this hybrid 
model by rationalizing iee resistanee estimates against other modeling teehniques, full 
power test results, and realistie performance of modern iee eutters. 
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II. BACKGROUND 


A. DEFINITIONS 

To more appropriately convey the subject material, all pertinent variables are 
hereby defined in Tables 1-5. All equations obtained from literature or developed are 
based on these definitions to aid consistency in presentation. Table 1 provides the unit 
vectors for two body-fixed coordinate frames. The icebreaker reference frame is given in 
Cartesian coordinates. The icebreaking surface frame is defined in two dimensions and is 
coupled to the icebreaker reference frame. The parameters and dynamics for the 
icebreaker are provided in Table 2. All parameters are typically known after the initial 
design stage. The icebreaker body has six degrees of freedom; three translational 
motions, and three rotational motions. This study emphasizes forward translation of the 
icebreaker only within this context, although follow on study should include additional 
motions. The parameters of sea ice including dimensional, material property, and material 
failure considerations are provided in Table 3. All physical parameters and variables for 
existing analytical ice resistance methods are included in Table 4. All parameters 
affecting the development and error of the hybrid method for average ice bending 
resistance are included in Table 5. An illustration of the reference frames, hull 
parameters, and motions are provided in Figure 1. 


3 



Table 1. Body-fixed eoordinate system eonventions. 


Reference Frame 

Coordinate System 

Unit Vector 

Definition 

Icebreaker Body-Fixed 

X 

longitudinal 


Y 

transverse 


Z 

vertical 

Hull-Ice Contact Surface 

n 

normal 


t 

tangential 


Table 2. Icebreaker parameters and motions variable definitions. 


Icebreaker 

Parameters 

Variable 

Name 


A 

displacement 


M 

mass 


Ps 

shaft power 

Dimensions 

L 

length at waterline 


B 

breadth 


T 

draft 

Bow Angles 

a 

waterline entrance 


3 

flare 


V 

stem 



effective ice contact 

Dynamics 

Variable 

Motion Affected 

Translation 

X 

surge 


y 

sway 


z 

heave 

Rotation 

0 

roll 



pitch 



yaw 
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Table 3. Sea iee parameters. 


Sea Ice 

Parameters 

Variable 

Name 

Dimensions 

d 

grain diameter 


H 

thickness 


Ic 

characteristic length 


1 

cusp length 


w 

cusp width 


R 

cusp radius 


Ro 

cusp radius reference 


n 

cusp radius angular position 

Properties 

E 

Young's modulus 


pice 

density 


V 

Poisson's coefficient 



friction coefficient 


va 

air concentration 


vb 

brine concentration 


Q 

temperature 

Material Failure 

ob 

bending stress 


ab,o 

uncorrected bending stress 


oc 

crushing stress 


TS 

shear stress 


ol 

principle stress 


a2 

principle stress 


do/dt 

stress rate 


6ol/5x2 

stress gradient 


C 

stress rate constant 


D 

stress gradient constant 


K 

stress gradient limit 
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Table 4. Physical parameters and resistance component variable definitions. 


Resistance Calculation 

Parameters 

Variable 

Name 


pw 

sea waterdensity 


Fp 

propulsion force 


Fv 

vertical ice loading force 


g 

gravitational acceleration 


n 

numberof cusps 


V 

velocity 


Vi 

initial velocity iterate 


Vi+1 

product velocity iterate 

Resistance 

Rb 

bending resistance 


Rc 

stem crushing resistance 


Rice 

ice resistance 


Row 

open water resistance 


Rr 

rotation resistance 


Rs 

submergence resistance 


Rsf 

skin friction resistance 


Rw 

wavemaking resistance 


Table 5. Hybrid method parameters for average ice bending resistance. 


Finite Element Bending Resistance 

Parameters 

Variable 

Name 


A 

FEM proportionality constant 


b 

error slope-intercept 


F 

distributed ice loading force 


m 

error slope 


N 

parameter power constant 


P 

substitute parameter 
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B. LITERATURE REVIEW 

I. Sea Ice Properties 

Accurate ice resistanee models require knowledge of sea iee material properties 
and failure mechanisms in the presenee of an ieebreaking ship. However, the material 
properties of sea iee are difficult to assign because many studies identify variations in iee 
samples based on age, ehemieal eomposition, temperature, seale, and testing method [9]. 
For example, in the Lindqvist 1989 study [13], nearly one order of magnitude variation in 
the flexural failure stress of iee was observed during in-situ testing. Therefore, the 
material properties obtained from literature are approximate quantities, and likely require 
refinement based on further speeification of these variables. 

Sea ice is formed from aeeumulated snow at the surfaee of seawater, and has a 
mierostrueture that differs eonsiderably from freshwater iee. Sea iee typieally assumes a 
hexagonal eolurnn strueture with anisotropie material properties. Unlike freshwater iee, 
whieh forms a uniform erystalline strueture, the grains in sea iee orient in several 
direetions, lending to the formation of an iee sheet with varying thiekness. Additionally, 
brine ehannels form between the grains, forming a heterogeneous solid solution within 
the iee sheet, and giving the iee viseoelastie properties. Repeated tests indieate either 
duetile or brittle material failure depending on the iee temperature, eomposition, and 
stress rate [9]. 
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a. Density 


-s 

The density of sea iee is roughly approximated as 910 kg/m [9]. 

b. Young’s Modulus 

Given the viscoelastic behavior of sea ice, Young’s Modulus varies from 
0.3 GPa-10 GPa for static loads, to 6 GPa-10 GPa for dynamic loads [9]. 

c. Poisson’s Ratio 

The Poisson’s ratio is approximately 0.295 [9], although it is also 
presented as 0.3 [9]. 

d. Flexural Strength 

(1) Temperature 

The flexural strength is related to the temperature of the ice for 

temperatures -0.4 C > > -18.5 C [9]. 

a-, =4.7-0.960-0.310' (2.1) 

(2) Composition 

The expression for flexural strength as a function of composition is 
shown below, where Va is the concentration of air, and Vb is the brine concentration [9]. 

( 2 - 2 ) 

(3) Stress Rate 

The flexural strength of the ice has no dependence on stress rate in 
the range 49.0 kPa/s < < 981 kPa/s. For flexural stress rates above 981 kPa/s, the flexural 
strength is obtained as a function of the stress rate [9]. 

^, dcT . 

^t=^b,o+C\n— (2.3) 

dt 


8 



(4) Testing Method 

Depending on either a small beam or ring tensile test for iee, the 
tensile stress of iee was measured at 0.785 MPa and 1.77 MPa, respeetively [9]. 

(5) Grain Diameter 

The flexural stress of iee is approximated as a funetion of the grain 
diameter, typieally between 1 em and 2 em [9]. 

1.11 


= 


d 


(2.4) 


Where the grain diameter, d is given in eentimeters, and the flexural stress, at is in MPa. 
e. Shear Strength 

The shear test estimates the shear failure stress at 2.06 MPa [9]. 


/. Compressive Strength 


Compressive failure of the iee varies by orientation of the iee erystal, 
having a maximum of 8.34 MPa for eores loaded on the prineiple axis, and a minimum of 
3.43 MPa for eores loaded perpendieular to this prineiple axis [9]. 


g. Stress Gradient 


In addition to these three failure meehanisms, sea iee fails at a stress 
gradient threshold [10]. 


(o-j - 0-2)" + o-j" + o-/- D 


ydX^j 


(2.5) 


This equation modifies the von-Mises plane stress failure eriterion to 
inelude the stress gradient term, the partial derivative of the prineiple stress with respeet 
to the perpendieular axis. 


2. Icebreaking Process 

Icebreakers induce failure of the ice by three observed mechanisms; crushing, 
shearing, and flexure. In most icebreaking operations, all three of these mechanisms are 
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present to some degree. Crushing and shear failure are localized at the hull-ice contact 
area, whereas flexural failure occurs at some distance from the icebreaker. The general 
pattern of damage generated by an icebreaker is a crushing zone at the stem of the 
icebreaker, with the formation of characteristic semi-elliptical cracks, known as cusps, 
from the centerline to the extreme breadth, caused by flexural failure. The cusps may also 
exhibit secondary crack formation from the contact surface radially along the length of 
the cusp [13]-[16]. 


Cusp Formation 



As a result of the material properties of sea ice, this pattern of cusp formation is 
dependent on icebreaker dynamics, as observed in laboratory and scale-model testing 
[15]. In certain conditions, the ice cusps occur predictably at specific locations along the 
bow. However, past research cannot conclude whether icebreaker dynamics are beneficial 
or detrimental to the overall ice resistance [15]. In addition, the interplay of failure 
mechanisms is dependent on the dimensions of the ice due to stress gradient limitations 
[10]. Upon failure of the ice, the icebreaker clears the ice underneath and to either side. In 
this stage, the icebreaker both rotates and submerges the broken pieces, providing some 
resistance to the ship. However, as the broken pieces are no longer providing structural 
resistance from the ice sheet, the thrust exceeds the resistance forces, and the icebreaker 
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regains momentum. To maintain transit speed through the iee, all momentum lost by the 
loading of the iee must be restored within a distanee equal to the length of the broken 
pieees in the direetion of motion, before the ieebreaker onee again eneounters the iee 
sheet [16]. 


3. Model Techniques 

Sinee icebreaking ships traditionally represent signifieant capital investment, 
modeling is conducted to allow iterative performance improvements to the intended 
design. Three types of modeling were developed, each with inherent advantages and 
disadvantages depending on the application. The three types of models are; scale models, 
analytical models, and numerical models. 

a. Scale Models 

Similar to tow tank testing, this modeling technique incorporates a layer of model 
ice floating on top of the water. The desired hull form is constructed, instrumented, and 
physically tested while breaking model ice. This technique is advantageous because 
resistance data is obtained directly, according to scale, and involves physical failure of 
the sea ice. The major caveat to this advantage is that the stress gradient, stress rate, and 
complex composition of sea ice disrupt the accuracy of the scale results. To more 
accurately simulate the full scale icebreaking process, the model ice is either “doped” by 
adding chemicals to affect the material properties, or is manufactured similarly to natural 
sea ice by allowing artificial snow to freeze on the tow tank surface [15]. Despite realistic 
interpretation of full scale icebreaking operations, the disadvantage to this method is that 
it is not practical at the initial design stages because of cost, time, and facilities 
availability. 


b. Analytical Models 

In contrast to scale modeling, analytical models [13], [14] are useful to 
calculate ice resistance at the initial design stages because of their simplicity. For ships 
following the general specifications of a modem icebreaker, existing analytical models 
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are the most praetieal means to determine iee resistanee. Sinee analytieal models are 
eonstrained by a eertain set of assumptions for a partieular method, a major disadvantage 
of a semi-empirieal analytieal model is that results are sealed with respeet to existing 
similar ships analysis, and may not provide reasonable performanee estimation for more 
radieal ship designs. Even with a standard hull form, in the ease of the iee eutting ship, 
eurrent analytieal methods are eonstrained to loeal iee failure, and eannot aeeount for the 
effeet of damage to the iee sheet at some distanee from the ieebreaker. However, one trait 
of analytieal methodology is the deeoupling iee resistanee into various constituents. For 
example; open water resistance, bending resistance, crushing resistance, rotational inertia 
resistance, and submergence resistance may be solved for separately, according to 
approximate system parameters, and scaled to match empirical data. Therefore, it is 
feasible that some modification of particular resistance constituents within the analytical 
method may allow inclusion of ice sheet damage to the overall ice resistance calculation. 

c. Numerical Models 

While finite element models have the potential to replicate the physical 
phenomena of breaking ice, these models encompass a wide range of programs with 
varying complexity and accuracy. Numerical methods have already implemented variable 
parameters and dynamics affecting icebreaking processes [17]-[19]. However, the 
inherent imperfections and nonlinearity of the solution of the numerical model requires 
validation with empirical data across the entire range of expected operations, adding to 
the cost of the technique. Although numerical models may provide insight into new 
radical icebreaker designs at marginal expense, they may also contain significant 
unforeseen errors that must be accounted for by other means. 

4. Cusp Geometry 

Several research efforts attempt to relate ice cusp dimensions to the ice 
loading because of its significant influence on ice resistance predictions. As a result, three 
types of ice cusp approximations were obtained from literature, with varying application 
and accuracy describing this form of damage realistically. In order of complexity, ice 
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cusps were analyzed according to one-dimensional, two-dimensional orthogonal, and 
two-dimension radial coordinate systems as shown in Figure 3. Within this speetrum of 
analysis, the ice is either uniformly loaded along the edge, or experiences point loading. 
The relative accuracy of each depends on the geometry of the hull. Since uniform loading 
is expeeted of bows with flat contaet surfaces, such as the Thyssen-Waas bow, the 
uniform loading analysis can reasonably deseribe the ieebreaking pattern. For other bow 
forms with some curvature, however, point loads are a more realistic assumption. In all of 
these cases, the ice is assumed as a homogeneous and isotropic sheet of uniform 
thickness with perfectly vertical edges on a uniform elastic foundation. 



a. One Dimensional 

The least complex form of ice cusp analysis assumes the ice sheet as a 
semi-infinite plate loaded along a single free edge. In this form, the ice is assumed to fail 
at a set distance from the free edge of the plate, proportional to the characteristic length 
of the plate, /c [13]. 
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This distance from the edge to the loeation of failure is referred to as the eusp length. Sinee 
this method fails to aeeount for the eusp width, two issues arise that affeet its applieation. 
First, it is diffieult to derive a realistie physieal assoeiation between the assumed distributed 
load required for eusp failure and the point load it represents. Seeond, the eusp length 
determined by this method is likewise restrieted to one-dimensional applieations. 
Therefore, sueh analysis is best suited to oases exhibiting: 1) distributed loading of the ioe 
sheet, and 2) the width of the distributed loading is muoh greater than the eusp length. For 
the majority of ioebreakers, however, this method is not reliable. 


b. Two Dimensional Orthogonal 

There are two variations of this method, oonstituting the progression from a 
distributed load to a point load assumption. This is aooomplished by assuming speoifio 
failure shapes aooording to presoribed geometries, in whioh the orthogonal eusp width is 
related to the eusp length. The eusp shapes are defined by either reetangular or triangular 
elements, as shown in the figure. In all eases, the eusp length is determined by the same 
method as the one-dimensional analysis. It is important to note that several analytieal 
ieebreaking resistanee methods eombine the reetangular and triangular elements, 
therefore, aeeepting a signifieant departure from both the natural ieebreaking proeess and 
analytieal deseription. 

(1) Reetangular Element 

This variation is similar to the one-dimensional method, exeept that the eusp 
width is related to the eusp length. This allows the frequeney of ice eusp formation to be 
determined, thus improving the aeeuraey of the total ioe load approximation, with 
marginal inorease in oaloulation. However, as in the one-dimensional ease, it is diffieult 
to assign an appropriate loading foroe that bears physieal oompatibility. 

(2) Triangular Element 

This variation relates the length of the eusp as sides of an isosoeles triangle to an 
arbitrarily determined base, the eusp width. This allows the eusp to be solved in the same 
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manner as the reetangular eusp, with exeeption of variable eross-seetion along the eusp 
length. However, sinee the free edge is infinitely small, the distributed foree aeting on 
this representative eusp may be approximated as a point load. In addition, these triangular 
eusps may be joined side by side to form a eusp resembling a semi-eirele experieneing a 
point load at the eenter. If an infinite number of triangular elements are used, the error 
between the eusp shape and triangular representation ean be effeetively eliminated. 
Therefore, this variation allows both versatile and somewhat realistie deseription of the 
eusp. However, the solution complexity limits it to numerical solvers [20]. 

c. Two Dimensional Radial 

The most complex, and likewise, accurate methods for determining the 
relationship between the ice loading force and semi-elliptical cusp dimensions 
appropriately analyze the system in a two-dimensional radial system. Similarly to the 
triangular cusps, the ice is loaded at a single point. Since the relationship obtained from 
literature has empirical basis, it is much more easily calculated upon determining the cusp 
variables. The caveat to this, however, is that a numerical solver is required to determine 
all variables [18]. 
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III. LINDQVIST METHOD 


A. REASON FOR COMPARISON 

The Lindqvist method for ealeulation of ice resistance is a semi-empirical 
analytical model that is referenced by several modern ice resistance methods as a 
standard for comparison [16], [21], [22], The Lindqvist method likewise serves as a 
benchmark for this particular study, although follow on research indicates that the 
method likely has some conservative error. Therefore, the hybrid ice resistance method is 
developed out of respect to the Lindqvist method, adopts components of this method to 
reduce error in comparison, and attempts to make improvements to experimental error. 

B. ASSUMPTIONS 

Several assumptions were made in the development of the Lindqvist method that 
affects the accuracy of the solution. These assumptions fall into two general categories, 
those that are in place to eliminate nonlinearity from the method, and those that provide 
rough approximation for poorly understood quantities. The empirical basis and linearity 
of the method allow it to be reasonably accurate, despite room for improvement within 
the individual assumptions. The assumptions of the Lindqvist Method are as follows: 

I. Linear Assumptions 

• Ice resistance increases linearly with velocity. 

• The vessel is motionless with exception of constant forward translational 
velocity. 

• The hull is approximated as a plane, developed from average hull angles, 
exerting normal and tangential force components on the ice at the 
waterline. 

• Upon contact, the ice fails instantaneously. 

• The failed ice pieces do not provide rotational inertia resistance back to 
the hull, nor are they affected by the viscosity of the water. 

• The ice does not exhibit viscoelastic, stress rate, or stress gradient 
properties. 
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2. Quantity Assumptions 

• Broken ice provides frictional resistance evenly distributed along 70% of 
the underside of the hull once submerged. 

• The resistance force is independent of the crushing failure stress. 

• The shear failure stress equals the bending failure stress. 

• Young’s modulus for ice is 2 GPa. 

• Poisson’s ratio for ice is 0.3. 

• Friction coefficient during hull-ice interaction is 0.1 for hulls with low- 
friction coatings and 0.16 otherwise. 

• Both the length and width of a failed ice cusp are equal to one third of the 
characteristic length of an edge-loaded semi-infinite plate with uniform 
elastic support. 

C. RESISTANCE COMPONENTS 

1. Stem Crushing 


The equation derived to approximate the stem crushing resistance is provided 
below. The equation is multiplied by an approximate proportional constant obtained by 
instrumenting the bow of the icebreaker with strain gauges, and taking transient 
measurements with reference to the ice bending failure stress and thickness. It is 
important to note that this resistance is not influenced by the dimensions of the 
icebreaker. 


Rc = 0.5 
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(3.1) 


where the effective ice contact angle C between the hull and the ice, is defined as a 
function of the waterline entrance and stem angles. 
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2. Bending 

The bending resistanee equation assumes the failure of a non-integer number of 
ice cusps with dimensions proportional to the characteristic length of the ice sheet, 
modeled as a semi-infinite plate supported on an elastic foundation. 

^ J //COSf ^ f. 1 
l^sinacosi^ cos^ 

3. Submergence 

The submergence resistance component provides the frictional contribution of the 
broken ice flows as a function of the approximated underwater surface area of the hull 
and submerged depth. This formulation incorporates the snow thickness resting on top of 
the ice because this is also submerged. This equation does not explicitly consider the 
influence of the added mass of water surrounding the ice chunks as they are submerged. 

T(^^]+f,k (3.4) 

\n + Al y 

where 
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D. EMPIRICAL VELOCITY RELATION 


As Lindqvist observed, measured ice resistance appears to increase linearly with 
velocity from some nonzero resistance value as expressed below. 


K,=(K+/t,) 1 + 1-4 


+ 3;. 1-I-9.4- 


The Lindqvist empirical velocity relation accounts for this observation by setting the 
zero-velocity resistance as the sum of the three resistance components, and increasing the 
value of each component linearly according to velocity, and scaling these values 
according to non-dimensional constant values. Thus, this velocity relation implies that 
two criteria must be met. First, assuming a linear relation between ice resistance and 
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velocity as observed, the three resistance components must equal the predicted zero- 
velocity slope-intercept of the empirical resistance data. Second, the contribution of each 
resistance component must be scaled appropriately with velocity, such that the sum of 
three resistance components equals the total resistance for a range of vessels. 

E. INTRINSIC VALIDATION 

The Lindqvist method was validated against both the calculated and empirical 
resistance data extrapolated from the original study. This comparison was necessary to 
ensure that the method was implemented properly into this study, since the Lindqvist 
resistance is typically presented as the total resistance, without specifying the 
contributions of the underlying resistance components. As desired, the resistance values 
calculated by Lindqvist were successfully reproduced for each vessel in the original 
study. Therefore, the individual resistance components analyzed by Lindqvist during the 
formation of his method were likewise reproduced, as shown in Figure 4. 



Figure 4. Relative influence of Lindqvist resistance components on the icebreaker 

Vladivostok with respect to velocity. 
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This analysis provided several insights into limitations of the method. Although 
the final empirical velocity relation achieves reasonably accurate calculation of the total 
ice resistance, there may be substantial errors in the individual resistance components. 
For example, the average resistance component values were obtained and compared with 
the total resistance versus vessel velocity. This indicated that although ice resistance 
appears to increase linearly with velocity, the contributions of the stem crushing and 
bending resistance components increase, while the submergence resistance component 
decreases. In addition, the relative contributions of each of these three resistance 
components were compared. Averaging from 0 m/s to 5 m/s, ice submergence provides 
approximately 40% to the total resistance, stem crushing provides 35%, and bending 
contributes 25%. Therefore, these percentages allow initial approximations of the ice 
resistance for icebreakers designed to eliminate any one component. 



Figure 5. Extrapolated raw data for the icebreaker Vladivostok comparing the original 
Lindqvist ice resistance calculation with a calculation using the method and 

the best fit line of the data. 
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Using the extrapolated data, and assuming a linear resistance velocity 
relationship, a best fit line was drawn from the raw data, as shown in Figure 5, since this 
was not provided in the original study. Overlaying the Lindqvist calculation over the raw 
data and best fit line revealed that although the Lindqvist data reasonably approximates 
the resistance, the zero-velocity resistance, nor the slope of the resistance line are 
properly matched by the method. This was expected from the assumptions of the method 
and analysis over a range of vessels, however, these errors may also stem from the 
absence of a pertinent resistance component, such as rotation of the ice. 

F. SENSITIVITY ANALYSIS 

Successful implementation of the Lindqvist method enabled sensitivity analysis to 
determine the most influential variables on total ice resistance with regard to the 
individual resistance components, located in Appendix A. Table 5 lists the parameters 
used for the sensitivity study of each ice resistance component. 

Table 6. Resistance component influential parameters 
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Lindqvist acknowledges that some portions of the method require improvement, 

therefore this analysis contains errors. Since Lindqvist also indicates that his method is 

more accurate for larger vessels, the largest vessel of the original study, Vladivostok, was 

selected to mitigate these errors. As expected, the sensitivity analysis determined that the 

ice resistance is highly dependent on the vessel gross dimensions, bow angles, and the ice 
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dimensions. However, the Lindqvist method proved highly sensitive to the bending stress 
of the iee and the density differential between sea water and sea iee. This dependenee on 
environmental variables requires that iee resistance calculations adequately describe the 
intended operating environment. Although the Lindqvist study described the bending 
stress for all cases, it did not describe the density differential. After several iterations, the 
best error obtained from adjusting the density differential was to within 10% of the 
original study. Therefore, the exact value of the density differential was not reproduced, 
but was close enough to proceed with analysis. 
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IV. HYBRID METHOD DEVELOPMENT 


The proposed hybrid method attempts to evaluate both eonventional and other 
ieebreaking strategies by aehieving synergy among all modeling methods. The ereation of 
this tool is motivated by the desire to perform objeetive ieebreaking resistanee analysis 
under a eommon framework, as shown in Figure 6. This framework is intended to 
provide the context to combine resistance components from multiple sources to 
determine the overall resistance. Comprehensive analysis of this method is not possible 
within the scope of this thesis because of its many variations across the range of all 
resistance approximation methods. Therefore, the proposed system requires further 
analysis, refinement of individual resistance components, and their relations. 

A. PROPOSED FRAMEWORK 

1. Resistance Components 

Research was conducted to determine all analytical resistance components from 
literature, listed below. 

• Bending 

• Rotation 

• Skin Friction 

• Stem Crushing 

• Submergence 

• Wave Making 

2. Icebreaking Functions 

Separation of resistance components was accomplished by first subdividing the 
icebreaking process into three necessary functions; breaking the ice, clearing the ice, and 
making headway. 


a. Icebreaking 

The stem crushing and bending resistance components are included within 


the icebreaking function because these resistance components are only applied as the 
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icebreaker comes into contact with the ice sheet although these resistance components are 
often summed as steady values, they are actually cyclical depending on the failure of the 
ice due to the failure of the ice extending some distanee ahead of the ieebreaker. The 
stem erushing eomponent provides more steady resistanee to the ieebreaker since the 
stem is the first part of the ieebreaker to contact the ice, and is relatively unaffeeted by 
cusp failure. The bending component eyclical frequency, on the other hand, is highly 
dependent on the dimensions of the ice cusps. 

b. Ice Clearing 

The clearing of broken ice involves rotation and translation of the broken 
chunks. This function is accomplished by the rotation and submergenee resistanee 
components existing in literature. The rotation ice resistance component is a function of 
both the rotational moment of inertia of the failed iee cusp and the icebreaker velocity. 
The submergence resistance component ineorporates both the translation of the broken 
ice under the water and the hull-ice contaet frietion, however, these may also be 
separated. 


c. Hydrodynamic Resistance 

Progressing relative to the iee sheet also involves making headway relative 
to the water beneath the ice. This relative motion ereates drag on the icebreaker body in 
the form of frietional and pressure drag. The skin frietion resistanee eomponent is direetly 
proportional to the surfaee area of the hull in eontaet with the water. The pressure 
differential generates waves through this medium at the water-air interfaee. This drag 
component is defined in literature as the wave making resistance. Since the skin friction 
and wave making resistance components are traditionally of great importance ship 
design, they are typically obtained as a coupled resistance from scale model tow ta nk 
testing, but may be aecounted for separately by a wide range of other methods. 
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Figure 6. Proposed icebreaker resistance framework. The comers of the triangle 

represent the three icebreaking functions. The Venn-diagram illustrates the 
icebreaking resistance components obtained from literature. The arrows 
indicate potential relationships between resistance components. 


B. ANALYTICAL VALIDITY 
I. Resistance Coupling 

The success of this method is governed by the validity of the assumption that the 
icebreaking process is comprised of discrete resistance components. This assumption is 
not particularly radical since the resistance components within the framework were all 
obtained from literature, and have been implanted into existing resistance models 
successfully [13], [16], [18]. Flowever, Figure 6 also identifies cases, both described in 
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literature and newly added, in which some resistance elements are coupled. In the broader 
sense, this illustrates the interdependence of the three icebreaking functions. The first 
such coupling was observed in previous research, that hydrodynamic forces resulting 
from the motion of the icebreaker contributed to failure of the ice by affecting the 
buoyant force supporting the ice sheet [16]. This has also been observed during high 
speed icebreaking by hovercraft [6]. The remaining three instances of coupling are 
speculated in this study. The first speculates that the bending and rotation resistances are 
linked by the geometry of the failed cusp. The second potential relationship links the 
submergence and hydrodynamic skin friction components by the finite wetted hull area. 
The third assumes that there is some connection between stem crushing and bending 
components based on the dimensions and form of the bow. 

2. Coupling Mitigation 

This lack of information provided the motivation to speculate which parameters 
influence this coupling such that analytical errors are reduced or expected. 

a. Velocity 

The wave making and skin friction resistance values become trivial as 
velocity approaches zero. Since these values also contribute to coupling of the bending 
and submergence components, respectively, it is reasonably assumed that analytical error 
is reduced as velocity approaches zero. The zero-velocity assumption does not have 
serious implications on icebreaker operations since it produces the operating limits. 

b. Hull Geometry 

The stem crushing, bending, and rotation resistance components are 
effectively negated as the hull-ice contact angle approaches zero. However, these values 
remain coupled. Since the only design parameter not shared by all three of these 
components is the vessel breadth, this value may be adjusted to reduce expected 
analytical error. The stem crushing resistance component is independent of vessel 
breadth. Therefore, the stem crushing is effectively de-coupled from the bending 
resistance in cases where the vessel is narrow, or has an extremely large breadth. Two 
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vessels that may provide further insight into this assumption at eaeh extreme are the E- 
Craft and the oblique ieebreaker eoneept, respeetively. 

C. RESISTANCE COMPONENT SUPERPOSITION 

Despite the validity of assuming truly deeoupled resistanee eomponents, the 
models presented in literature approximate ieebreaking resistanee as the sum of either 
average or diserete values over the ieebreaking eyele. If a steady state resistanee value is 
desired at a speeified veloeity, the resistanee eomponents obtained from either means 
may be used interehangeably aeeording to the work-energy prineiple. 



In addition, these values may be interehanged to obtain a rough resistanee estimate for 
nonlinear systems. 

I. Conceptual Test 

To test this eoneept, the Lindqvist ealeulation was performed for the ieebreaker, 
Vladivostok. In the first test, it was assumed that the average value of all resistanee 
eomponents were present through the entire duration of the ieebreaking eyele. The 
steady-state veloeity of the ieebreaker was obtained via applieation of the work-energy 
prineiple. In the seeond test, the bending resistanee eomponent was assumed to oeeur 
only upon eontaet with the iee sheet until failure of the iee eusp, with all eusps failing 
simultaneously. To further test this eoneept, the seeond test was also modified sueh that 
that individual eusps were assumed to fail at equal intervals. 

a. Averaged Ice Resistance 

Upon applying the work-energy method presented in Equation 4.1 to the 
transient numerieal model, the ieebreaker aehieved steady-state veloeity as displayed in 
Eigure 7. This method ineludes the mass of the ieebreaker, and therefore, more aeeurately 
represents the veloeity profile of an ieebreaker. In this averaged iee resistanee ease, all 
resistanee values are assumed to be present at every veloeity iterate. In addition, the iee 
resistanee inereases proportionally with veloeity. Sinee the hydrodynamie resistanee was not 
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included in this model, the inerease in ice resistanee causes the eonvergence of the velocity. 
In this particular case, the steady state veloeity eonverged to approximately 0.325 m/s. 


Vladivostok: Continuous Cusp Failure Velocity 
Profile at Rjce = 600 kN 

0.35 



Figure 7. Continuous cusp failure transient velocity profde. 


b. Cyclical Bending Resistance 

The steady-state veloeity profde for simultaneous and single cusp failure 
was obtained using Equation 4.1, appearing as a saw-tooth pattern as expected from 
literature [13]. Figure 8 illustrates the contrast between assuming simultaneous and 
individual cusp failure. This model differs from the average ice resistance model because 
it has the added dependence on the absolute distance traveled, along with veloeity 
dependence. The distance traveled is significant in this model because it determines the 
locations of ice loading prior to cusp failure over suecessive icebreaking cyeles. The saw¬ 
tooth veloeity profile indieates nearly instantaneous momentum change of the ieebreaker. 
When all cusps are assumed to fail simultaneously, the model may predict a lower initial 
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velocity limit, below which, the icebreaker is unable to proceed through the ice. Although 
this highlights the motivation for backing and ramming icebreaking operations, where the 
icebreaker cannot proceed continuously, this assumption assumes over predicts the 
propulsive force required to transit the ice for a given initial velocity. To provide a more 
realistic minimum velocity at the start of the icebreaking cycle prior to the necessitation 
of backing and ramming, the velocity profile resulting from the failure of a single ice 
cusp was maintained. This is also shown in Figure 8 as a saw-tooth velocity profile 
having significantly less cyclical amplitude than that predicted with the simultaneous 
cusp failure assumption. 


Vladivostok: Failure of Single v. All Cusps Velocity Profile at 

= 600 kN 
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Figure 8. Cyclical cusp failure transient velocity profile. 


2. Influence on Total Resistance 

The steady-state velocity was converged upon for a resistance range of 534 kN, as 
expected for zero-velocity by the Lindqvist equation, to 1400 kN. As shown in Figure 9, 
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the continuous cusp failure assumption exactly overlaid the Lindqvist calculation. The 
cyclical cusp failure assumption, on the other hand, was unable to match neither the 
slope, nor the slope-intercept of the Lindqvist calculation. 


Vladivostok: Continuous v. Cyclical Numerical Ice 
Resistance Calculation 



Velocity [m/s] 


. Lindqvist Calculation 

♦ Continuous 
▲ Cyclical 


Figure 9. Comparison of continuous and cyclical cusp failure with icebreaker 

Vladivostok raw data and Lindqvist resistance predictions versus velocity. 


The cyclical cusp failure calculation error was not expected, because this should 
have perfectly matched the Lindqvist calculation line. Flowever, the average error 
between this numerical calculation to the values expected from the Lindqvist calculation 
over the resistance range was less than three percent, as shown in Figure 10. The likely 
cause of this error is from processing limitations. To assist this argument, the steady state 
velocity was obtained with a MATLAB® code included in Appendix B. Inspection of the 
command window output indicated that computational errors developed from calculating 
the square root term of the velocity iterates. 
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Figure 10. Comparison of continuous and cyclical cusp failure resistance-velocity 
relationship against icebreaker Vladivostok raw data. 


D. CONCEPTUAL TEST DETERMINATION 

The hybrid methodology is most easily evaluated when coupling of the resistance 
components is assumed negligible. By assuming zero velocity, thereby negating the 
rotation and hydrodynamic resistance components, the hybrid method takes a form 
resembling the Lindqvist method. To test the concept of multi-sourced superposition, a 
resistance component was obtained from a second model, and compared with the 
Lindqvist component. The bending resistance component was chosen for comparison 
since the Lindqvist approximation has a non-zero value at zero-velocity and since it is the 
only resistance component that is not localized to the hull-ice contact zone. Therefore, 
testing the resistance component has the three-fold purpose as a proof of concept for the 
hybrid method, cross-examination of the Lindqvist bending resistance, and to provide the 
ability to quantify the effects of damage to the ice analytically. 
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E. 


FINITE ELEMENT MODEL 


Given the zero-veloeity assumption to reduee potential analytieal errors, the 
bending resistanee was obtained from a finite element interpretation of the ieebreaker 
providing a statie distributed load to the iee sheet at the bow of the ieebreaker, as shown 
in Figure 11. 



Figure 11. Illustration of iee loading, eusp length, and transformed eusp length 
determination in statie finite-element iee bending model. 

Similar to existing analytieal methods, the iee sheet is assumed homogeneous and 
isotropie, with no spatial or time dependent properties, and supported on a uniform elastie 
foundation with stiffness equal to the produet of sea water density and gravitational 
aeeeleration. The finite element model was generated using ANSYS Statie Struetural 
toolbox. Given the limitations of finite element modeling, it was desirable to seleet 
dimensions for the iee sheet sueh that the iee behaved as a semi-infinite plate eonsistent 
with prior analytieal methods and that the iee had a suffieient number of nodes elements 
through the thiekness of the iee sheet to eliminate numerieal errors as the iee sheet was 
deformed. Assuming the ieebreaking proeess is symmetrieal about the eenterline of the 
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icebreaker, only one half of the ice sheet was modeled to increase caleulation speed. The 
dimensions of this half sheet of ice were arbitrarily seleeted as 100 m long by 35 m wide. 

lee was removed from one edge along the length of the ice sheet inwards to 
approximate the area oceupied by the half-breadth of the icebreaker in question. In addition, 
this cut assumed a triangular bow form, aceording to the average waterline angle, to roughly 
approximate the influenee of this parameter on the bending resistanee. The ice thickness was 
varied through a realistie range of 0.1 m to 2.0 m. The iee sheet was bounded on the aft, 
forward, and outboard edges by the fixed boundary eondition. The symmetric boundary 
condition was used along the centerline edge of the ice sheet. A uniform elastic support was 
applied on the underside of the ice sheet. A distributed load was applied at the hull contact 
area to roughly estimate the loading on the ice. To improve accuracy of the solution, mesh 
refinement was performed in the vieinity of the hull eontact area using three-dimensional 
tetrahedral elements. The solution output was the bending failure stress gradient versus the 
distributed load aeting on the iee, ineluded as Figure 12. 



Figure 12. FEM statie structural bending failure of the iee modeled in ANSYS. 


The distance of ice failure ahead of the ieebreaker bow was measured at multiple 
points with respeet to the icebreaker breadth. This distanee was obtained from multiple 
ice thicknesses and a one order magnitude range for the iee bending stress, from 0.1 MPa 
to 1 MPa. 
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F. ANALYSIS 

As both the loading force and distance of broken ice ahead of the icebreaker are 
obtained from the finite element model, this is simply related to the existing analytical 
formulations for the bending resistance. Two obstacles develop, however, since analytical 
ice bending failure models neither represent the force as a distributed load, nor do they 
represent the distance ice is broken ahead of the icebreaker. The latter of the two is 
accounted for since it is the transformed cusp length. The implications of distributed 
loading, on the other hand, required further insight. As expected, without accounting for 
this inconsistency, the relation between the finite element and Lindqvist predictions for 
the bending resistance was not clear across the range of ice thicknesses. Critical 
observation of the Lindqvist method revealed that it assumed both cusp length and width 
to be equal independent of ice thickness. However, other literary sources indicated that 
the length and width were in fact dependent on ice thickness. Since Lindqvist 
approximates the bending resistance as the sum of point loads, spaced apart by the cusp 
width, this total resistance component is proportional to a cusp width dependent 
approximation. Figure 13 illustrates the load distribution and cusp failure geometry 
assumed as a single icebreaking cycle. 
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Observed 




Figure 13. Illustration of ice loading and cusp failure geometry from real-world 

observation, the Lindqvist 2-D orthogonal ice cusp failure assumption, and 
the static structural finite element ice failure assumption. 


Therefore, the finite element approximation for the average bending resistance 
force is proportional to the Lindqvist value. 


F 


Fn 


T_ 

= A 

1 


FEM 
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(4.2) 


where A is the proportionality constant. 

The average bending force predictions obtained by finite element method and 
Lindqvist were compared according to parameter variation through a reasonable range of 
values. In this comparison, the Young’s modulus, bending failure stress, and thickness of 
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the ice, as well as the breadth of the vessel, are clearly defined in the Lindqvist equation, 
whereas the influence of these parameters is not clear within the finite element result. The 
range of parameter variation is provided in Table 7. 


Table 7. Parameter variation range 


Parameter 

Min 

Max 

E 

0.2 Gpa 

20 Gpa 

Ob 

0.1 Mpa 

1.0 Mpa 

H 

0.1m 

2.0 m 

B 

10 m 

30 m 


Due to the enhanced ability of the finite element model to approximate the cusp 
failure geometry, some error was expected between the finite element and Lindqvist 
average force predictions. The error was linearized in all cases by multiplying the error 
by the variable parameter to some power. This linearization defined the proportionality 
constant. A, between the finite element and Lindqvist average force predictions. The 
average bending force obtained from the finite element result was then directly applied as 
a substitute to the Lindqvist value. 


^ m-p + b'' 


(4.3) 


Although this form of analysis has limitations, it generally had good agreement for the 
given parameters. The slope and slope-intercept of the linearized relative error was 
obtained and provided in Table 8. 


Table 8. Parameter modification variables for proportionality constant, A. 


P 

N 

m 

b 

E 

0.095 

4.00E-08 

-602.74 

ob 

0.185 

-0.006 

-41.42 

H 

-0.015 

40.406 

-70.215 

B 

0.9 

-69.964 

968.6 


All graphs in this section are presented together in Appendix C as a quick 
reference to this statement. 
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1. Young’s Modulus 

The Young’s Modulus was varied from 0.2 GPa to 20 GPa to aeeommodate the 
wide range of values presented in literature, sinee this variable is expeeted to ehange 
signifieantly in the expected operating environment. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 14, showing nonlinear 
dependence on parameter variation. 


Uncorrected Relative Error; 
Ice Young's Modulus Variation 



-120 - 

Young's Modulus [Pa] 


Figure 14. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice Young’s modulus. 
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b. Linearized Relative Error 

The data was linearized by multiplying the error by the Young’s modulus 
to the power of 0.095, produeing a best fit line with an R value of 0.998, as shown in 
Figure 15. 


Linearized Relative Error; 
Ice Young's Modulus Variation 



-700 

Young's Modulus [Pa] 


y = 4E-08x-602.74 
= 0.998 


- FEM Error 

-Linear (FEM Error) 


Figure 15. Linearized relative error of the finite element average bending foree v. 

Lindqvist when varying the iee Young’s modulus. 


c. Corrected Average Bending Force 

The best fit error predietion line was then used to eorreet the finite element 
predietion to the Lindqvist value aeeording to Equation 4.3. The eorreeted average 
bending foree obtained by finite element method is eompared with the Lindqvist 
prediction in Figure 16. The data sets are consistent, as expected from this analysis. 
However, the slope-intercept of the finite element corrected value appears to be 
significantly greater than that predicted by Lindqvist. Since the average bending force 
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curves are nearly vertical as the Young’s modulus approaches zero, this error was 
expected by numerical processing limitations. 



Figure 16. Corrected finite element average bending force v. Lindqvist when 

varying the ice Young’s modulus. 


d. Corrected Relative Error 

To finalize the parameter variation of Young’s modulus, the relative error 
of the corrected finite element prediction for the average bending force was compared 
with the Lindqvist value, as shown in Figure 17. As stated before, significant errors were 
expected as the Young’s modulus approached zero. However, if this region is neglected, 
the finite element obtained force nearly matehes the Lindqvist prediction with a 
maximum error of 11.3% and a minimum of 1.5%. 
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Corrected Relative Error; 
Young's Modulus Variation 


70 



Figure 17. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice Young’s modulus. 

2. Bending Failure Stress 

The error between the Lindqvist and finite element average bending force 
predictions was obtained for the bending stress range of 0.1 MPa to 1 MPa to be 
consistent with the range of values obtained during in-situ measurement. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 18, showing nonlinear 
dependence on parameter variation. 


42 











Figure 18. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice bending failure stress. 


b. Linearized Relative Error 

The data was linearized by multiplying the error by the bending failure 
stress to the power of 0.185, producing a best fit line with reasonable agreement to the 
errors presented, as shown in Figure 19. 
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Figure 19. Linearized relative error of the finite element average bending force v. 

Lindqvist when varying the ice bending failure stress. 


c. Corrected Average Bending Force 

The best fit error prediction line was then used to correct the finite element 
prediction to the Lindqvist value according to Equation 4.3. The corrected average 
bending force obtained by finite element method is compared with the Lindqvist 
prediction in Figure 20. The data sets are consistent, as expected from this analysis. 
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Figure 20. Corrected finite element average bending force v. Lindqvist 
when varying the bending ice failure stress. 


d. Corrected Relative Error 

To finalize the parameter variation of the bending failure stress, the 
relative error of the corrected finite element prediction for the average bending force was 
compared with the Lindqvist value, as shown in Figure 21. Although the error curve is 
nonlinear, the maximum and minimum errors are 2.3% and 0.3%, respectively. 
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Corrected Relative Error; 

Ice Bending Failure Stress Variation 


5 



Figure 21. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice bending failure stress. 

3. Ice Thickness 

The error between the Lindqvist and finite element average bending force 
predictions was obtained for ice thickness ranging from 0.1 m to 2.0 m. This variable is 
expected to change significantly in the expected operating environment. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 22, showing limited if any 
dependence on parameter variation. 
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Uncorrected Relative Error; 
Ice Thickness Variation 



-70 

Ice Thickness [m] 

Figure 22. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice thickness. 

b. Linearized Relative Error 

Although the prediction appears to fit well without correction, the data 
was linearized by multiplying the error by the bending failure stress to the power of - 
0.015, producing a best fit line which mitigated residual errors, as shown in Figure 23. 
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Linearized Relative Error; 
Ice Thickness Variation 


20 



y = 40.406x-70.215 
R* = 0.9989 


- FEM Error 

-Linear (FEM Error) 


Figure 23. Linearized relative error of the finite element average bending force v. 

Lindqvist when varying the ice thickness. 


c. Corrected Average Bending Force 

The best fit error prediction line was then used to correct the finite element 
prediction to the Lindqvist value according to Equation 4.3. The corrected average 
bending force obtained by finite element method is compared with the Lindqvist 
prediction in Figure 24. The data sets are consistent, as expected from this analysis. 
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Corrected Average Bending Force; 
Ice Thickness Variation 



0 0.5 1 1.5 2 

Ice Thickness [m] 


Lindqvist 

— — Corrected FEM 


Figure 24. Corrected finite element average bending force v. Lindqvist when varying 

the ice thic kn ess. 


d. Corrected Relative Error 

To finalize the parameter variation of the ice thickness, the relative error 
of the corrected finite element prediction for the average bending force was compared 
with the Lindqvist value, as shown in Figure 25. Despite the curve presented in the error 
plot, the maximum and minimum errors were 1.0% and 0.3%, respectively. 
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Corrected Relative Error; 
Ice Thickness Variation 



Figure 25. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice thickness. 

4. Vessel Breadth 

The error between the Lindqvist and finite element average bending force 
predictions was obtained for vessel breadth in the range of 10 m to 30 m, consistent with 
most icebreaking vessels. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 26, showing nonlinear 
dependence on parameter variation. 
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Figure 26. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the vessel breadth. 


b. Linearized Relative Error 

The data was linearized by multiplying the error by the bending failure 
stress to the power of 0.9, producing a best fit line with reasonable agreement to the 
errors presented, as shown in Figure 27. 
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Figure 27. Linearized relative error of the finite element average bending force v. 

Lindqvist when varying the vessel breadth. 


c. Corrected Average Bending Force 

The best fit error prediction line was then used to correct the finite element 
prediction to the Lindqvist value according to Equation 4.3. The corrected average 
bending force obtained by finite element method is compared with the Lindqvist 
prediction in Figure 28. The data sets are consistent, as expected from this analysis. 
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Figure 28. Corrected finite element average bending force v. Lindqvist when varying 

the vessel breadth. 


d. Corrected Relative Error 

To finalize the parameter variation of vessel breadth, the relative error of 
the corrected finite element prediction for the average bending force was 
compared with the Lindqvist value, as shown in Figure 29. Despite the curve 
presented in the error plot, the maximum error was 0.11% and the minimum error 
was 0.09%. 
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Corrected Relative Error; 
Vessel Breadth Variation 

0.25 -- 



Figure 29. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice vessel breadth 
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V. DAMAGED ICE SHEET RESISTANCE PREDICTION 


A. TEST SPECIFICATION 

Following the development of the hybrid bending resistance calculation, the finite 
element model of the ice sheet was altered to represent specific damage patterns prior to 
interaction with the icebreaker hull. This modified ice sheet was then loaded with the 
icebreaker hull, and this load was substituted for the analytical bending resistance 
component in the Lindqvist calculation. As the substitution errors are relatively small for 
the range of values selected in this initial study, the variation in total ice resistance 
resulting from the damage to the ice sheet contains minimal error within the limitations of 
the static structural finite element model. Although there are infinite ice damaging 
strategies that may be implemented, this study only entertains the concept of linear cuts 
travelling parallel to the icebreaker at three locations along the breadth of the vessel; at 
the vessel centerline, at Im outboard of the extreme breadth, and 2m outboard of the 
extreme breadth. The cuts were defined as either partial penetration through 50% of the 
ice sheet, or as full penetration cuts and were assumed to occur at some distance well 
ahead of the icebreaker. In total, 8 specific damage patterns were tested. The tests are 
identified in Table 9. 


Table 9. Ice damage test case identification 


Test Case 

Location 

Depth 

1 

Side [Im] 

50% 

2 

Side [Im] 

100% 

3 

Side [2m] 

50% 

4 

Side [2m] 

100% 

5 

Center 

50% 

6 

Center 

100% 

7 

Side [Im] and Center 

50% 

8 

Side [Im] and Center 

100% 
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B. RESULTS 


The ice was damaged according to the patterns defined in the test specification, 
and this was accounted for as an averaged bending force. The finite element bending 
failure stress contours are included in this section to illustrate the location of damage and 
the expected failure regions. The obtained bending resistance component for damaged ice 
relative to undamaged ice is included in Table 10. 


Table 10. Relative bending resistance of undamaged ice versus damaged ice 

by test case. 


Test Case 

Location 

Depth 

Rb [%] 

1 

Side [Im] 

50% 

132.40 

2 

Side [Im] 

100% 

11.16 

3 

Side [2m] 

50% 

119.54 

4 

Side [2m] 

100% 

147.52 

5 

Center 

50% 

159.23 

6 

Center 

100% 

24.63 

7 

Side [Im] and Center 

50% 

68.17 

8 

Side [Im] and Center 

100% 

6.00 


1. Side Cut [Im] 

This cut involves a cut Im outboard of the extreme breadth of the icebreaker, 
running parallel to the vessel’s longitudinal dimension in the finite element model. 

a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thickness. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model; however, the stress contour is visible at the cut location. The finite element 
solution to this damage pattern for the bending failure stress contour is included as Figure 
30. 
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Figure 30. Bending failure stress contour for a partial side cut located Im outboard 

of the extreme breadth of the icebreaker. 


b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cut. The ice failure region differs considerably from the undamaged ice model, indicating 
stress concentrations at both the forward most extent of the cut and at the icebreaker 
shoulder. The finite element solution to this damage pattern for the bending failure stress 
contour is included as Figure 31. 
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Figure 31. Bending failure stress contour for a full penetration side cut located Im 

outboard of the extreme breadth of the icebreaker. 
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2. Side Cut [2m] 

This cut involves a cut 2m outboard of the extreme breadth of the icebreaker, 
running parallel to the vessel’s longitudinal dimension in the finite element model. 


a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thickness. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model, however, the stress contour is visible at the cut location. The finite element 
solution to this damage pattern for the bending failure stress contour is included as Figure 
32. 
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Figure 32. Bending failure stress contour for a partial side cut located 2m outboard 

of the extreme breadth of the icebreaker. 
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b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cut. The ice failure region differs considerably from the undamaged ice model, indicating 
stress concentrations at both the forward most extent of the cut and at the icebreaker 
shoulder. The finite element solution to this damage pattern for the bending failure stress 
contour is included as Figure 33. 
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Figure 33. Bending failure stress contour for a full penetration side cut located 2m 

outboard of the extreme breadth of the icebreaker. 

3. Center Cut 

This cut involves a cut at the centerline of the icebreaker, and running parallel to 
the vessel’s longitudinal dimension in the finite element model. 

a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thickness. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model, however, a high stress region extends forward of the stem of the icebreaker at the 
cut location. The finite element solution to this damage pattern for the bending failure 
stress contour is included as Figure 34. 
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Figure 34. Bending failure stress contour for a partial penetration centerline cut located 

at the stem of the icebreaker. 


b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cut. The ice directly ahead of the icebreaker failed similarly to the undamaged ice model; 
however stress concentrations occurred at both the forward most extent of the cut and at 
the icebreaker shoulder. The finite element solution to this damage pattern for the 
bending failure stress contour is included as Figure 35. 
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Figure 35. Bending failure stress contour for a full penetration centerline cut located 

at the stem of the icebreaker. 
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4. 


Side and Center Cut 


This cut involves both a cut Im outboard of the extreme breadth of the icebreaker 
and at its centerline, running parallel to the vessel’s longitudinal dimension in the finite 
element model. 


a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thickness. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model; however, high stress regions develop at the shoulders and directly forward of the 
stem of the icebreaker at the cut locations. The finite element solution to this damage 
pattern for the bending failure stress contour is included as Figure 36. 



Figure 36. Bending failure stress contour for partial penetration side and centerline 

cuts, located at the extreme breadth and stem of the icebreaker, respectively. 


b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cuts. The ice failure region differs considerably from the undamaged ice model; 
indicating stress concentrations at both the forward most extent of the cut and at the 
icebreaker shoulder. The finite element solution to this damage pattern for the bending 
failure stress contour is included as Figure 37. 
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Figure 37. Bending failure stress contour for full penetration side and centerline cuts, 
located at the extreme breadth and stem of the icebreaker, respectively. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

1. Hybrid Method 

A hybrid method was developed sueh that iee resistanee data from multiple 
sourees may be eross-refereneed or eombined under eommon framework to better 
deseribe ieebreaking in the analytieal eontext. Within this framework, a numerieal 
approximation for the bending iee resistanee eomponent was obtained via finite element 
model and related to existing analytieal model with empirieal basis in order to better 
understand the underlying physies as well as to improve the aeeuraey of existing iee 
resistanee models. Upon aehieving low error for given predefined influential parameters, 
the finite element iee sheets were damaged aeeording to speeified damage patterns, and 
this damage was aeeounted for as a deerease in the total iee resistanee by using the 
Lindqvist analytieal model. 

2. Damaged Ice Sheet Resistance 

The finite element model of the iee sheet was altered aeeording to eight damage 
seenarios to relate any ehange in the bending stress failure eontour to the bending 
resistanee eomponent. The method presented in this study indieates that the bending 
resistanee is signifieantly redueed when the iee is damaged by either full penetration outs, 
or outs extending ahead of the ioebreaker stem. The reduotion of bending resistanee for 
the oenterline out is partioularly interesting beoause this is also the looation of the 
breadth-independent stem crushing component. It may be possible to reduce the stem 
crushing component as well as the bending resistance for centerline cuts in the ice. Since 
the stem crushing and bending resistance components are likely coupled, this study was 
unable to account for any reduction in the stem crushing resistance while solving for the 
bending resistance. The present finite element static structural model could not account 
for stress gradient or stress rate failure criteria for the ice sheet because these values are 
not well represented in literature. 
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B. FUTURE WORK 


To proceed with this work, both the time and size dependent failure properties of 
sea ice must be obtained. Without this information, it was impossible to numerically 
model the transient failure of ice and test against empirical data. In addition, further tests 
must be performed to quantify resistanee eomponent eoupling in the hybrid model 
framework so that analytical methods may be further refined. If, for instance, the 
eoupling between the stem erushing and bending resistance eomponents are resolved, the 
effeets of ice damage on the stem erushing component may be pursued. 
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APPENDIX A. LINDQVIST SENSITIVITY ANALYSIS 


Each parameter in the Lindqvidst resistance formula was varied across a 
reasonable range of values, with respect to the vessels used by Lindqvist to formulate the 
method. The figures below illustrate the effect of this variation on the individual 
resistance components. 
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APPENDIX B. NUMERICAL ICEBREAKING VELOCITY 

MATLAB CODE 


close all 
clear all 
clc 

% Icebreaker Vladivostok 
% Parameters 
L = 112; 

B = 23.5; 

T = 9.5; 

alpha =18; % waterline entrance angle 
gamma =24; % stem angle 

zeta = 56; % effective ice contact angle 
mu = 0.16; 

H = 0.58; 

Hsnow = 0.28; 

E = 2E9; 
nu = 0.3; 
sigmab = 500; 
sigmac = 3430; 
rhow = 1026; 
deltarho = 109; 

m = 5000; 
g = 9.81; 

Fpmin = 0; 

Fpmax = 1400; 

Fv = 84.1; 

Rc = 67.25344; 

Rb = 105.383; 

Rs = 360.9507; 

Ic = (E*H"3/(12*(l-nu"2)*rhow*g))"0.25; 

1 = lc/3; 

lx = 1/sind(alpha); 

Ixn = 10; 

1x1 = Ix/lxn; 

a = Fv/ (sigmac*10''3*B) ; 
ax = a/sind(alpha); 
axn = 10; 
axi = ax/axn; 

s = ax+lx; 
sn = axn+lxn; 

n = B/(l*sind(alpha)); 
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Fbmax = Rb*(ax+lx)/(2*ax); 

for k = Fpmin:Fpmax 
Fp = k; 

Fprop(k-Fpmin+1,1) = Fp; 
i = 1; 

Vn = 1.5; 

Vs = 0; 

Vsf = 1; 

while Vsf-Vs > 0.05 
for j = 1:axn 


Vs = 

Vn; 




Fb = 

Fbmax/axn* j; 




Rci 

= Rc*(1+1.4*Vn/sqrt 

(g 

*H) ) 

r 

Rbi 

= Fb*(1+1.4*Vn/sqrt 

(g 

*H) ) 

f 

Rsi 

= Rs*(1+9.4*Vn/sqrt 

(g 

*L) ) 

r 

Vi (j 

, i) = Vn; 




Vo = 

Vi(j,i); 




Vx = 

2*axi/m*(Fp-Rci-Rb 

i-: 

Rsi) 

+Vo^2; 

Vi (j 

+ 1, i) = sqrt(Vx); 




Vn = 

Vi(j+1,i); 




V( j + 

(i-1)*(sn+1)) = Vn; 




j = 

axn+1:sn+1 




Vi (j 

, i) = Vn; 




Fb = 

0; 




Rci 

= Rc*(1+1.4*Vn/sqrt 

(g 

*H) ) 

f 

Rbi 

= Fb*(1+1.4*Vn/sqrt 

(g 

*H) ) 

f 

Rsi 

= Rs*(1+9.4*Vn/sqrt 

(g 

*L) ) 

r 

Vo = 

Vi(j,i); 




Vx = 

2*lxi/m*(Fp-Rci-Rb 

i-: 

Rsi) 

+Vo^2; 

Vi (j 

+ 1, i) = sqrt(Vx); 




Vn = 

Vi(j+1,i); 




V( j + 

(i-1)*(sn+1)) = Vn; 





end 

Vsf = Vn; 
i = i+1; 

end 

Vss(k-Fpmin+1,1) = k; 

Vss(k-Fpmin+1,2) = mean(Vi(:,i-1)); 

end 

plot(V) 

%Lindqvist = [0 534;4.251 1400]; 
%plot(Vss(:,2),Vss(:,l)) 

%hold on 

%plot (Lindqvist (■.,!), Lindqvist ( : , 2) ) 
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APPENDIX C. HYBRID METHOD DEVELOPMENT ERROR 


Linearized Relative Error 
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